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ABSTRACT 



A plane, heated jet of a constant property fluid exiting 
into a plane, straight channel of the same fluid is analyzed. 
The effect of buoyancy on the velocity and temperature fields 
is investigated. Two jet orientations are considered; a 
transverse jet, which enters the stream perpendicular to the 
main flow, and a longitudinal jet, which enters the stream 
parallel to the main flow. 

Solutions for the temperature and velocity fields for low 
Reynolds number flow are obtained using a finite difference 
scheme. Results are presented for three different values of 
the buoyant force, holding other variables constant. For an 
isothermal jet, results for the velocity field are presented 
at a higher Reynolds number. 
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NOMENCLATURE 



List of 


Symbols 


a 


Thermal diffusivity 


c 


Specific heat 


Ci 


Coefficients of interpolating polynomial 


D 

•N. 


Channel depth 


D j 


Width of jet 


Dj* 


Non-dimensional width of jet 

/ 


F 


Unknown general function 


g 


Acceleration of gravity 


Gr 


Grashof number based on channel depth 


h 


Coefficient of heat transfer 


H 


Grid spacing in vertical direction 


J (A , B) 


Jacobian operator 


k 


Thermal conductivity 


K i 


Coefficients of finite difference approximation 


M 


Square matrix 


P 


Pressure 


P (x,y) 


Interpolating polynomial 


Pr 


Prandtl number 


q 


Heat flux 


Re 


Reynolds number based on channel depth 


T 


Temperature 


T o 


Undisturbed stream temperature 


u 


Velocity in horizontal direction 
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u 


Non-dimensional velocity in horizontal direction 


u. 

D 


Non-dimensional jet velocity 


U o 


Average undisturbed stream velocity 


U s 


Surface velocity 


V 


Velocity in vertical direction 


V 


Non-dimensional velocity in vertical direction 


V i (x , y) 


Variable terns of interpolating polynomial 


X 


Space coordinate in horizontal direction 


X 


Non-dimensional space coordinate in horizontal 
direction 


y 


Space coordinate in vertical direction 


y 


Non-dimensional space coordinate in vertical 
direction 


a 


Ratio of grid spacing in horizontal direction 
to spacing in vertical direction 


6 


Coefficient of thermal expansion 


e 


Non-dimensional temperature 


e £ 


th 

Solution for 0 at 1 iteration 


u 


Absolute viscosity 


V 


Kinematic viscosity 


p 


Density 


>p 


Stream function 


Ip* 


Non-dimensional stream function 


lp* z 


Solution for at Jt*** 1 iteration 


v 2 


Harmonic operator 


v 4 


Biharmonic operator 


< > 


Row vector 


{ } 


Column vector 


C ] 


Square matrix 


Cr 


Order of 
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List of Subscripts 



a 

D 

i 

j 

jL 

ju 

L 
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s 

T 

u 



Ambient 

Downstream boundary 
Incoming stream 
Jet 

Lower edge of jet 
Upper edge of jet 
Longitudinal jet 
Maximum 
Surface 

Transverse jet 
Upstream boundary 
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I. 



INTRODUCTION 



Every commercial power generation plant releases some 
waste heat to the environment. Water from a nearby stream, 
lake, ocean is normally used as a reservoir for this 
rejected heat. As the demand for power increases, more and 
more heat will be discharged into these bodies of water. 

The advent of nuclear power plants, which have a lower ther- 
mal efficiency than conventional fossil-fueled plants will 
tend to increase the amount of waste heat per unit power 
output. With the emergence of these trends, concern about 
the effect of these thermal discharges has begun to grow. 

The utilization of a body of water for cooling purposes 
is normally carried out by diverting a certain mass flow 
from the body of water, passing it through a heat exchanger, 
and returning it at an increased temperature. Hence, the 
temperature of the body of water, at least in the local 
sense, is increased. 

From the biological viewpoint, there is no doubt that 
raising the temperature of water can be harmful to some forms 
of aquatic life. For example, for some fish a seasonal 
temperature rise of as little as 2°C can be fatal. Aside 
from such obvious damage, raising the temperature of water 
decreases its capability to retain dissolved oxygen, increases 
the metabolic rate of fish, and affects their reproductive 
habits and ability to resist disease. In addition, an 
increase in environmental temperature is beneficial to some 
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potentially harmful forms of plant life and bacteria [ 1 ] . 
Although the immediate effect of such changes may not be 
obvious, the normal equilibrium of the ecological system 
has been disturbed, and the possibility of long range 
damage does exist. 

In order to ensure that temperature disturbances are 
within safe limits, the precise nature of the temperature 
disturbance must be known. The quantitative value of the 
increase in temperature above the undisturbed temperature 
of the given body of water as a function of position and 
time would be of significant value in analyzing the effect 
of a thermal discharge. Several studies related to this 
problem have been carried out. 

The overall effect of a thermal discharge on a given 
body of water is one method of approaching the problem. A 
method for predicting the equilibrium temperature of a 
cooling lake was described by Sefchovich [2], A study 
relating plant parameters to the average upstream tempera- 
ture in a stream or river was developed by Nahavandi and 
Campisi [3], 

An experimental investigation of the temperature profiles 
in the vicinity of a cooling plant located on Yankee River 
was presented by Merriam [4], Approximate analytical models 
were compared with various experimental temperature profiles 
downstream of power plants situated on rivers by Polk, 
Benedict, and Lahey [5], 

A possible model of a thermal discharge is a heated jet 
exiting into a given environment. For the case of thermal 
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discharge into a lake or ocean this model takes the form of 
a heated jet exiting into a quiescent fluid. An analytical 
solution for the two dimensional temperature and velocity 
profiles for the case of a heated, vertical , laminar jet 
exiting into an infinite environment was developed by Brand 
and Lahey [6], Additional work on a heated jet exiting 
into a quiescent fluid is now being conducted at Oregon 
State University [7]. 

The case of a heated discharge into a stream or river 
may be modeled by a heated jet exiting into a moving stream. 
An experimental investigation of a buoyant, turbulent, 
circular air jet exiting into a cross wind was conducted by 
Ramsey and Goldstein [8], An analytical study of the 
surface spread of a submerged, buoyant incompressible jet 
was made by Tulin and Schwartz [9]. 

The problem undertaken in this study was that of a 
heated jet in a stream. This particular problem can be 
approached in several ways. The general problem is a tran- 
sient problem in which the temperature and velocity fields 
vary in three dimensions. Factors which must be considered 
are diurnal temperature and heat flux variations, tidal 
effects, seasonal temperature and flow variations, atmos- 
pheric conditions, stream meander path, bottom contour 
variations, fluid property variations, the possible turbulent 
nature of the effluent, and temperature stratification 
effects. Because of the complexity of this problem, several 
simplifications were made. 
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The transient conditions were assumed to be "slow" 



compared to the response time of the stream. Therefore, 
at any instant in time the stream could be considered to 
be at a quasi-equilibrium state. The time dependent terms 
were therefore not considered. The fluid properties were 
assumed to be constant, for simplification, and the axis of 
the river was assumed to be straight. 

It was decided to attempt to isolate the effects of. 
buoyancy on the temperature and flow patterns. Also, to 
put the problem in a more tractable form, it was decided to 
consider variations in two dimensions only. Assuming that 
one of the dimensions would be the downstream dimension, 
there are two possible choices: a plane parallel to the 

river surface or a vertical plane passing through the axis 
of the river. Since the buoyancy forces act in the vertical 
plane, the vertical plane was chosen. In order to isolate 
buoyant effects, the bottom of the channel was assumed to be 
straight and parallel to the river surface. Also, for sim- 
plification, the river and the jet were assumed to be 
laminar in nature. The desired solutions were the tempera- 
ture and velocity profiles as a function of position and 
significant non-dimensional parameters. 
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II. ANALYSIS 



A. COORDINATE SYSTEM 

The problem under consideration in this work is flow 
in a channel into which a known quantity of fluid is injected. 
The injected fluid is at a higher temperature than the fluid 
in the channel. Since it is desired to study the effects of 
buoyancy, variations in the vertical plane only are consid- 
ered. The axis of the river is assumed to be straight. 
Rectangular cartesian coordinates are used, as shown in 
Figure 1. 

Two different orientations of the injected mass with 
respect to the main flow are considered. The longitudinal 
jet is injected parallel to the direction of main flow. 

The transverse jet is injected perpendicular to the direction 
of main flow. For the longitudinal jet, the origin of the 
coordinate system is located at the channel bottom, immedi- 
ately below the jet location. For the transverse jet, the 
origin of the coordinate system is located far enough 
upstream of the jet so that the incoming temperature and 
velocity profiles are not disturbed (see Figure 1). 

B. GOVERNING EQUATIONS 

Applying the Boussinesq approximation, which assumes 
that all density variations other than those giving rise to 
the buoyant forces may be neglected, and neglecting viscous 
dissipation, the steady two dimensional forms of the conti- 
nuity, Navier-Stokes , and energy equations are: 
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where T Q 


is the 


incoming 


stream 


temperature . 





By cross-differentiating with respect to x and y, and 
by subtracting equation (3) from equation (2) , the pressure 
p is eliminated. The resulting equation is: 



9x 2 9x9y 9x 2 9y 9y J 9x' > 9x9y z 9x 



+ v — 

9x9y 9y^ 



+ aV A_£v_ 3T 



(5) 



The continuity equation may be eliminated by defining 
a stream function ^ , such that 



u = 


9 ^ 

9 y 


, V = 


94> 

9X 






Only two unknown quantities remain 


. to be 


determined , 


^ and T. The 


two 


remaining 


equations 


are 
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Equation (6) is the momentum equation expressed as the 
transport of vorticity, where the vorticity is defined as 

V 2^ . 

The equations will be non-dimensionalized by introducing 
the following variables: 



U = 



V = 



v_ 

V. 



X 



x 

D 



Y = 



Z 

D 



e 




ip * 






UqD 



where U Q is the average undisturbed velocity, D is the 
channel depth, and Tj is the jet temperature. 

By substituting these into equations (6) and (7) , 



dip* 

dY 



2 

9 V ip* 
[ — ] 



3^* dV 2 ip* 

— [ 

9 X 9 Y 




Gr 3 e 
Re 2 9 X 



( 8 ) 



ill i! _ ill 91 = _±_ v 2 e 

3 Y 9X 9X 9 Y 



RePr 



(9) 



where Re is the Reynolds number based on the channel depth, 

Pr is the Prandtl number of the fluid, and Gr is the Grashof 
number based on the channel depth and the temperature differ- 
ence between the jet temperature and the incoming stream 
temperature . 

C . BOUNDARY CONDITIONS 

Examination of the governing equations reveals that 
twelve boundary conditions are necessary to completely 
specify the problem. The necessary conditions are four 
conditions for i() on x, four for on y, two for T on x and 
two for T on y. These conditions will be specified on the 
channel bottom, the river surface, upstream and far down- 
stream of the jet. The boundary conditions for \p are first 
expressed in terms of velocity, and then related to the 

stream function from the definition of \p , i.e., = u, 

3y 

3 ip _ _v 

9x 

At the upstream boundary, the undisturbed incoming 
horizontal velocity profile is assumed to be a 20% truncated 
parabola, which is a common assumption in stream flow [10]. 

A 20% truncated parabolic flow profile has the maximum 
velocity located two tenths of the total depth below the 
free surface (see Figure 2) . For the longitudinal jet, 
this profile is modified by assuming that the stream 
velocity is equal to the jet velocity across the jet width 
(see Figure 3) . For the transverse jet, the upstream 
boundary is located far enough upstream so that the velocity 
profile is equal to the undisturbed velocity profile. .The 
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Figure 2 - 20% truncated parabolic velocity profile. 



y 




Figure 3 - Velocity profile at upstream boundary 
for longitudinal jet. 
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vertical component of the stream velocity is assumed to be 
zero. For the longitudinal jet, the stream temperature is 
assumed to be equal to T 0 except across the jet width, 
where the temperature is equal to the jet temperature, T j . 
For the transverse jet, the upstream boundary is located 
far enough upstream so that the incoming temperature 
profile is also undisturbed. 

On the channel bottom, the no slip boundary condition 
requires that u = 0. The vertical component of velocity 
is also equal to zero, except at the transverse jet. In 
that case, v is equal to zero along the river bottom except 
across the jet width, where v is equal to the jet velocity. 
The river bottom is assumed to be adiabatic, yielding a 
temperature gradient boundary condition. 

On the channel' surface, it is assumed that there are 
no waves; i.e., v = 0. By assuming that the pressure is 
constant, and that the river surface is level, the appli- 
cation of the Bernoulli equation along the channel surface 
leads to the result that the surface velocity U s is 

constant. Since U is constant, and the value of U c at 

s “ 

the upstream boundary is known, the surface velocity is 
equal to the undisturbed surface velocity. At the surface, 
the conductive, convective, and evaporative nodes of heat 
transfer to the atmosphere are lumped into an overall heat 
transfer coefficient. It is assumed that radiative heat 
fluxes can be modeled by a known constant heat flux q. 
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A simple heat balance shows that the heat conducted into 
the surface plus the radiative heat flux into the surface 
equals the heat convected out. 

The velocity profile at the downstream boundary is 
assumed to be parabolic in form such that the total mass 
flow is equal to the undisturbed mass flow plus the mass 
flow of the jet. The vertical component of velocity is 
assumed to be equal to zero. It is also assumed that all 
heat added by the jet is convected out of the channel 
through the surface. Hence, the temperature profile 
returns to its original constant value T Q . The downstream 
boundary is situated far enough downstream so that the 
above conditions are met within a given tolerance. 

The mathematical forms of these conditions are: 
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where the L subscript refers to the longitudinal jet, the 
T subscript refers to the transverse jet, the i subscript 
refers to the incoming velocity profile, the D subscript 
refers to the velocity profile at the downstream boundary, 
yj L is the y-location of the lower edge of the jet, yjy 
is the y-location of the upper edge of the jet, Xj L and 
Xjy are defined similarly for the x-location of the jet, 
and T is the ambient air temperature. 

Since the problem is cast in terms of the stream 
function , it is desired to reduce the boundary conditions 
to a more convenient form. At the undisturbed upstream 
boundary, the 20% truncated parabola gives the following 
form for the velocity profile. 

u = |1 = ^[8^-5 (I) 2 ] 

ay 3 D D 

Integrating from 0 to D, it was found that the average 
velocity U c is given by 




Now, integrating u (y) from zero to y, assuming 

<J>(x, 0) = 0 

U~D „ _ 

* = — [12 (I) 2 - 5 (I) 3 ] 

r U 7 D D 

This is the upstream boundary condition for 4 > for the 
transverse jet. For the longitudinal jet, if it is assumed 
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that the jet diameter is small with respect to the river 
depth, the expression for ip becomes: 

ip = <f> u (y) °£y < YjL 

= My> + u j D j ( Y D ~ j j ~ L ,) YjL < y < y ju 

= 4» u (y) yj U < yi D 

where Dj is the jet width, yj L is the location of the lower 
edge of the jet, and is the location of the upper edge 

of the jet. 

For the longitudinal jet, the boundary condition at 
the channel bottom may be satisfied by assuming is con- 
stant. Hence, \p is chosen to be zero on the channel 
bottom for this case. 

For the transverse jet, the boundary condition for v 
may be integrated, giving 



i|> T = 0 for x<x jL 

x~x • 

*T = " U j D j ( “~ ) x jL <x<x jU 

+T = -u j D j x>x jU 

The boundary condition for v at the surface may also 
be satisfied by choosing ij> = constant. The value of the 
constant is equal to the value of ip at the surface as 
given by the upstream boundary condition. 

The downstream boundary condition for t|» may be found by 
integrating the downstream velocity profile, and evaluating 
constants from the surface and bottom values for ^ . 
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After performing the indicated operations, and non- 
dimensionalizing , the final form of the boundary conditions 
becomes 



X=0 



iKp*(Y) = ^ U *(Y) = 1 (12Y^ - 5Y J ) 



> 0 * (Y) 
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3 Grp 
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where the L and T subscripts have the same meaning as 
defined previously. 
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III. METHOD OF SOLUTION 



A. GOVERNING EQUATIONS 

The governing equations for the problem under considera- 
tion are: 

a ifj* a v^* a ijj* ^v 2 ^* 

3 Y 3X 3X 3 Y 

3 if)* 36 dip* 36 1 

— = V^e (11) 

3 Y 3X 3X 3 Y RePr 

In order to solve this problem using a finite differ- 
ence method, a finite difference operator for each 
differential operator in the governing equation is required. 
The method used to determine these operators is basically 
the same as that used by Dias [ll]. The operators developed 
in this study, however, are applicable to rectangular grid 
networks as well as to square grid networks. 

The procedure for determining the finite difference 
operators is to first define an interpolating polynomial 
P (x, y) having n undetermined coefficients. This poly- 
nomial can be used as an approximate to any general function 
F. (F can be interpreted as a stream function, for example.) 
By requiring that P (x, y) be equal to F at each of n 
points, n equations involving the n undetermined coeffi- 
cients are generated. Hence, the coefficients of P (x, y) 
can be determined in terms of the value of F at each node. 



= L.,4, 



,* - 



Re 



Gr 36 
Re 2 3 X 



(10) 
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By applying the desired differential operator to P (x, y) , 
an approximation for the operator in terms of the value of 
F at each node is obtained. 

To develop the biharmonic operator, a polynomial must 
be chosen such that the fourth partial derivatives with 
respect to both x and y are non-zero. The grid network 
through which P (x, y) is to be passed is shown in Figure 
4. Since there are 25 points in this grid, 25 unknown 
coefficients are required in the polynomial P (x, y) . 

Hence, an appropriate choice for P (x, y) is 

P(x,y) = C^y^x^ + C 2 y^x^ + Cgy^x^ + C 4 yx^ + C 5 X^ 

+ Cgy^x 2 + C^y^x^ + Cgy 2 x 2 + Cgyx 3 +Cio x3 
+C^^y^x^ +C^2Y^ x ^ +c 13Y 2x3 + C'14Y x ^ +C 15 x ^ 

+ c 1 6 Y 4x +Ci7Y 3x +Ci8 Y 2x +Ci9Y x +C 2 o x 

+ c 21Y 4 + c 22Y 3 +c 23^ 2 +c 24Y • +C 25 (12) 

This equation may be written in matrix form as: 

P(x,y) = <V i (x,y)> {C i } i = 1,2,. ..,25 (13) 

where (x, y) is x^y^ corresponding to Cj_ in equation (12) . 

This polynomial is now required to pass through the 
value of Fj at each node j, where Fj is the value of the 
function F at node j. At each node j, the following 
relation is obtained: 

P( x j/Yj) = <V i (x j ,y j )> (Ci) = Fj i = 1,2,. ..,25 

(14) 
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Figure 4-5x5 rectangular grid. 



By applying this relation at each of the 25 grid points. 



the following set of equations results: 
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(15) 

j = 1,2,..., 25 



Or, in matrix notation: 



{Fj} = [m] { C j } 



j = 1,2,...,25 



(16) 



Solving this matrix equation for Cj , 



{Cj} = [M] -1 {F j } 



j = 1,2, ... ,25 



The coefficients Cj of P (x, y) are now known, and 
P (x, y) can be used as an approximation for F. The 
biharmonic operator may now be applied to the polynomial 
P (x, y) . 

V 4 (P (x,y) ) = <V 4 ( V i (x,y)> (C^) 

= V 4 < (Vi (x,y) ) > (C i ) i = 1,2,. ..,25 

Since an operator for the biharmonic operator for the 
central node is desired, this expression is evaluated at 

x 13 '^13 * 

V 4 (P (x ,y ) ) = <V 4 (V i (x,y))> [M] _1 {F. } i=l,2,...,25 

1 (17) 

= <K ± > (F i > i=l , 2 , . . . , 25 (18) 
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where <K^> is a row vector of constants obtained by 
multiplying the matrix expressions appearing in equation 
(17) . The result is that are the coefficients of a 
finite difference approximation for the biharmonic operator 
applied at the central node of the grid in Figure 4. 

In order to ensure geometric consistency, this same 
interpolating polynomial was used to develop the operators 
appearing on the left hand side of equation (10) . The 
desired operator is applied to the polynomial, and evalu- 
ated at the central node. The results for a = 5 are 
presented in Figures 5, 6, and 7, where a is the ratio of 
the spacing in the x-direction to the spacing in the 
y-direction. 

A finite difference approximation for the harmonic 
operator, which is required in equation (11) is available 
in the literature [12]. Consistent approximation for the 
other operators appearing in this equation are also 
available. These operators are presented in Figures 8 
and 9 . 

B. BOUNDARY CONDITIONS IN FINITE DIFFERENCE FORM 

Assuming there is no jet present, the boundary condi- 
tions for \p* are all of the form 

ip* = f(t) 




where t is the space coordinate tangent to the boundary, 
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rence approximation for V operator for uniform 
grid (a=5) , applicable at central node. 
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rectangular grid (a=5) , applicable at central node. 
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Common Factor 




Truncation Error = ^(H 2 +o 2 h 2 ) 



Figure 8 - Finite difference approximation for 

V 2 operator applicable at central node 
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Figure 9 - Finite difference approximation for 

9/3X or 3/3Y applicable at central node 
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n is the coordinate normal to the boundary, and f (t) and 
g (t) are the values specified at each particular boundary. 
The first condition on \p* is satisfied in the finite 
difference scheme by evaluating f (t) at each node on the 
given boundary. These values are then used in the finite 
difference computations. 

To satisfy the boundary condition for the normal deriva- 
tive of rp , g (t) is evaluated at each node on the boundary. 

A finite difference approximation for the normal derivative 
is then required. An offset three point operator is used 
because the truncation error is consistent with the trunca- 
tion error of the governing equation. The finite difference 
operators used at the boundaries are presented in Figures 
10 and 11. 

The temperature boundary conditions assume various 
forms. At the upstream and downstream boundaries, the value 
of the temperature at the boundary is specified. At nodes 
located on these boundaries, the known value of the temper- 
ature is used in the finite difference computations. On 
the channel bottom, the derivative of the temperature is 
specified. At nodes on this boundary, the offset three 
point operator utilized previously is applied. At the 
channel surface, the sum of the temperature and its deriva- 
tive is specified. At nodes on the channel surface, a 
mixed boundary condition exists. In this case, the value 
of the temperature at the node located on the surface is 
added to the three point operator shown in Figure 11. 
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For most computations, the jet center is located between 
two nodes. In addition, the jet diameter is assumed to be 
less than the grid spacing in all cases. The effective 
diameter of the jet is equal to the grid spacing on the 
boundary on which the jet is located. This is a result of 
the fact that any length smaller than the grid spacing 
cannot be accurately modeled in a finite difference network. 
This can best be illustrated by example. Figure 12 depicts 
two possible jet configurations, both having the same volume 
flow rate. Also shown are the stream function profiles, 
superimposed on a grid network. In both cases, the jet 
width is less than the grid spacing. It can be readily 
seen that the values of the stream function at the two 
nodes on either side of the jet are the same, regardless 
of the jet diameter. The value of the stream functin is 
changed by the amount of the volume flow rate of the jet. 

If the jet width is smaller than the grid spacing, the 
effect of the jet on the finite difference grid is a fixed 
flow rate between two grid points. The effective diameter 
of the jet is then equal to the grid spacing, and the 
effective jet velocity is the flow rate divided by the grid 
spacing . 

By modeling the jet in this manner, the edges of the 
jet fall on nodes of the finite difference grid. In a 
continuous system, there is a step change in both the 
temperature profile and the velocity profile at these 
points. This step change cannot be modeled precisely in a 
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Figure 12 - Velocity profiles and resulting stream 
function profiles for two jets located 
between nodes 1 and 2. 
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finite difference scheme. By specifying the temperature 
at these two nodes in the finite difference grid, the 
temperature profile is as shown in Figure 13. It can be 
seen that the total area under the temperature profile 
curve is equal to the area under a step change of tempera- 
ture 2 x T s . The energy input of the profile shown in 
Figure 13 is equal to that of the step change of height 
2 x T s . This temperature may be defined as the effective 
jet temperature. 

If the three point finite difference operator shown in 

Figure 9 is used to calculate the velocity profile based 

on the relation V = 3 ^*. , it can be shown that the effec- 

3Y 

tive jet velocity defined earlier is consistent with the 
definition of the effective jet temperature. 

The above examples are directly applicable to the 
transverse jet. By analogy, the results may easily be 
extended to the longitudinal jet. 

For the results used to construct the convergence plot, 
the jet center is located at a node. For this orientation, 
the effective jet diameter is twice the grid spacing. 

C. NUMERICAL SOLUTION SCHEME 

In the numerical solution of this problem, the governing 
equations are rearranged 

v4^* = _ 11 + [ill 9 ( 72 4>*) _ iivV)_]Re (18) 

Re 3X 3 Y 3X ,3X 3Y 

v 2 e = RePr [ 3 ^ l 9 . - 3 ^ 3 9 ] ( 19 ) 

3 Y ax 3X 3 Y 
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Figure 13 - Temperature profile over finite difference 
grid in the vicinity of the transverse jet. 
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These equations may be expressed in an alternate form if 



we recognize that the Jacobian operator is defined as 



J (A , B) 



9A 9 B 9 A 9 B 



9X 9 Y 9 Y 9X 



Equations (1) and (2) may then be rewritten as: 




Gr 90 t/i* n 2 ... 
nr v? " J(^*,v ip* )Re 



( 20 ) 




RePr J(4<*,0) 



( 21 ) 



In the solution of these equations, an iterative scheme 
based on the uncoupled linear homogeneous forms of these 
equations is used. Once this linear system of equations is 
solved, the right hand side of equations (20) and (21) may be 
calculated . 

The linear equations are then resolved, using the calcu- 
lated values of the right hand side of the equations as known 
inhomogeneities. This procedure is repeated until the desired 
convergence criteria is satisfied. 

The uncoupled linear homogeneous forms of equations (20) 
and (21) are: 



The method of solving these equations with the boundary 
conditions for the transverse will be illustrated using the 
grid in Figure 14 as an example. 




0 



0 
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The boundary conditions on \ji* specify the value of ip* on 
all boundaries. Therefore, the value of ip* at nodes located 
on the boundaries is known. The value of ip* at each interior 
node is unknown. In order to solve for the value of ip* at 
interior nodes, one equation for each node is required. 

For all ip* values not adjacent to any boundary (nodes 15, 
16, 21, and 22 in example), an equation is obtained by applying 
the governing equation at the node. The finite difference 
operator shown in Figure 5 is used. 

For all ip* values located at nodes adjacent to the upstream 
boundary (nodes 8-11 in example) , an equation is obtained by 
applying the boundary condition I4 - = 0 at the adjacent 
boundary node (nodes 2-5 in example) . The finite difference 
operator shown in Figure 10 is used. 

For ip* values at nodes adjacent to the downstream boundary 
(nodes 26-29 in example) , an equation is obtained by applying 
the boundary condition |4— =0 at the adjacent boundary 
node (nodes 32-35 in example) . The finite difference operator 
shown in Figure 11 is used. 

For values at nodes adjacent to the channel bottom, but 

not adjacent to either the upstream or downstream boundary 

(nodes 12, 20 in example), an equation is obtained by applying 
the boundary condition -r-^— =0 at the adjacent boundary node. 

d x 

The finite difference operator shown in Figure 10 is used. 

For values at nodes adjacent to the channel surface, but 

not adjacent to either the upstream or downstream boundary, 

(nodes 17, 23 in example), an equation is obtained by applying 
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9 ] h * 9 

the boundary condition = =- at the adjacent boundary 

node (nodes 13, 19 in example). The finite difference operator 
shown in Figure 11 is used. 

Since the value of ip* is known at nodes located on the 
boundary, whenever a finite difference operator involves the 
value of tp* on a boundary, the known value of ip* is used. 

The resulting system of equations is a set of algebraic equa- 
tions for each unknown. This system of equations is solved 
by Gaussian elimination. 

The boundary conditions for 0 specify the value of 0 at 
the upstream and downstream boundaries, and at the two nodes 
on either side of the jet. However, the value of the tempera- 
ture on the channel surface is unknown, and the value of the 
temperature on the channel bottom, except for the nodes on 
either side of the jet, is also unknown. Again, one equation 
for each unknown is required. 

For each 0 unknown not located on a boundary (nodes 8-11, 
14-17, 20-23, 26-29 in example), the equation V^0 =0 is 
applied at the node at which the unknown is located. The 
finite difference operator shown in Figure 8 is used. 

For 0 unknowns located at nodes on the channel bottom 
(nodes 19-25 in example) , an equation is obtained by applying 
the boundary condition t-tt =0 at the node. The finite 

a X 

difference operator shown in Figure 10 is used. 

For 0 unknowns located at the channel surface, (nodes 
12, 18, 24, 30 in example), an equation is obtained by 
applying the boundary conditions ^+^2-0 = 0 at the node. 
The finite difference operator shown in Figure 11 is used. 
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The result is an equation for each 6 unknown. This system 
of equations is solved using Gaussian elimination methods. 

The method of solving equations (20) and (21) is presented 

in the form of a flow chart in Figure 15. The superscripts, 

o th 

e.g ., ip* refer to the solutions at every node on the Z 

iteration in a particular loop. The quantity R s (i|i* m ,0 n ) is 

defined as the right hand side of equation (20) based on the 

known values of ij;* 111 and 0 n . The quantity R Q is the right hand 

side of equation (21) defined similarly. 

The convergence tests for ip* and 0 are: 



. i-1 . i 

}P ~_JP 



T-T 



< 0.001 



e 1 ' 1 - e 1 



T-1 



< 0.001 



It is assumed that when the above criteria are satisfied 
at every node, the solution has converged. 

It should be pointed out that the governing equation is 
not applied at every unknown node. The non-linear terms 
calculated from the values of functions from the previous 
iteration affect the known value in the governing equation only 
at nodes at which the governing equation is applied. However, 
the solution for every node is affected, because the equations 
are solved simultaneously by Gaussian elimination techniques. 

In the calculation of the right hand sides of equations (20) 
and (21) , two different techniques were used. The first 
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Figure 15 - Flow chart for iterative scheme. 
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technique utilized the operators shown in Figure 6, 7, and 9, 
which were derived earlier. 

A second method , which provided improved numerical 
stability, was also used. This method is described in [13]-. 
This is simply an alternate way of calculating the Jacobian, 
which appears both in the right hand side of equation (20) 
and equation (21) . 

Using this method, the Jacobian Operator is expressed in 
three alternate forms: 



J (A , B) 



3 A 


3 B 3 A 


3 B 






3X 


3 Y 3 Y 


3X 






3 


(A |f) - 


3 


(A 


3 B 


3X 


3 Y 


3X 


3 


<*•£> - 


3 


(B 


3 A. 


3 Y 


3X 


3 Y 



When the quantities are grouped in this manner, applying 

the three point central difference operator shown in Figure 9 

leads to three different numerical methods of calculating 

J(A,B). By calculating J(A,B) by each of these methods, and 

using the average of the three values, some forms of numerical 

instability are eliminated. 

This method may be applied directly to equation (20) . 

However, in order to apply it to equation (21) , the value of 
2 

V ip* must be calculated at each node. This is done by applying 

the finite difference operator shown in Figure 8. The method 

2 

described above can then be applied to calculate J(i|»*,V t|>*) . 

The computer program used to obatin the results which are 
presented is included in Appendix B. 
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IV. RESULTS 



A. PRESENTATION OF RESULTS 

In the results, velocity data are presented in the form 
of streamlines in the channel. Temperature data are presented 
in the form of isotherms in the channel. Both the streamlines 
and the isotherms are presented in non-dimensional form, and 
are plotted with respect to non-dimensional values X and Y. 
These graphs are based on linear interpolation between calcu- 
lated values of ip* and 9 at nodes in the grid. It should be 
noted that in the graphical plots for both ip* and 0, the scale 
in the Y-direction is twice that in the X-direction. The grid 
used for all calculations consisted of 7 nodes in the Y- 
direction and 19 in the X-direction. 

For all results presented, the surface boundary condition 
on the temperature is ^ + 6 = 0 . Additionally, the 

ratio of the flow rate of the jet to the flow rate of the 
river is 0.2 in all cases. This value was chosen as a repre- 
sentative value for an actual power plant, and was felt to be 
small enough so that the surface boundary condition on the 
vertical velocity (i.e., no waves) remains valid. 

In Figures 16-21, the streamlines and isotherms for a 
transverse jet of width 5D/6 are presented. The Reynolds 
number and Prandtl for all these figures are Re = 0.3 and 
Pr = 6.0. Results for three different values of the ratio 
of the Grashof number to the Reynolds number are presented; 
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Gr/Re = 600, Gr/Re = 0, and Gr/Re = -300. The case of a 
negative value of Gr/Re implies that the temperature of the 
jet is less than the temperature of the stream. Hence, the 
buoyant force acts in the opposite direction, i.e., downward. 

In Figures 22-27, the streamlines and isotherms for a 
longitudinal jet of width of D/6 located between Y = 1/3 
and Y = 1/2 are presented. For all these figures. Re = 0.3 
and Pr = 6.0. The three values of Gr/Re presented are 
Gr/Re = 600, Gr/Re = 0, and Gr/Re = -300. 

In Figures 28-29, the streamlines for both a transverse 
jet of width 5D/6 and a longitudinal jet of width D/6 
between Y = 1/3 and Y = 1/2 are presented. The values of 
the non-dimensional parameters for these figures Re = 40, 

Pr = 6.0, and Gr/Re =0. No isotherms are presented because 
the jet is assumed to be at the same temperature as the 
river. 

In Figure 30 and Table 1, the surface temperature as a 
function of X is presented. The temperatures for the longi- 
tudinal jet are presented in tabular form becasue differences 
in the temperatures for the three different values of Gr/Re 
are difficult to distinguish graphically. For all cases in 
Figure 30 and Table 1 the Reynolds number and the Prandtl 
number are Re = 0.3 and Pr = 6.0. The value of Gr/Re for 
each case is indicated in the results. 

B. DISCUSSION OF RESULTS 

The effect of the buoyant force on the flow patterns and 
isotherms can be seen by examining the graphical results in 
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Figure 30 - Surface temperature vs. X for a transverse jet 
of width 5D/6, Re = 0.3, Pr = 6.0. 
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TABLE I 



Surface Temperature ( 0 ) 

o 

Longitudinal Jet 



X 


6 s 

Gr/Re = 600 


6 

s 

Gr/Re = 0 


6 s 

Gr/Re = 


0.0 


o 

• 

o 


0.0 


o 

• 

o 


0.833 


0.1083 


0.1077 


0.1074 


11. 667 


0.0907 


0.0901 


0.0898 


2'. 500 


0.0632 


0.0635 


0.0636 


3.333 


0.0431 


0.0437 


0.0440 


4.167 


0.0294 


0.0300 


0.0304 


5. 000 


0.0200 


0.0300 


0.0304 


5:. 833 


0.0137 


0.0141 


0.0144 


6.667 


0.0094 


0.0097 


0.0099 


7.500 


0.0064 


0.0066 


0.0067 


8.333 


0.0044 


0.0046 


0.0047 


9.167 


0.0030 


0.0031 


0.0032 


10.000 


0.0021 


0.0022 


0.0022 



-300 
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Figures IS— 21 ami 2ZH 2277.. Thee case of: Gf/Ree== PC implies the 



buoyant effect: is nagflliigj.hi'iBr: iir. the. momentUmr.equation;. . In 
this case,, ncr buoyant: farces aedbss on thee flow/ , andcthesmomentum 
equation iLs independent: off the; temper alureeequationr . A posi- 
tive value of Gr/Kte implies? that: the - temperature eof ithe jet is 
greater than T .. Iir. this case, the buoyant-- force eaets in the 
positive F-d r re ctiarr. . A 7 -. negatives values of:' Gf /Reeimplies that 

the temperature of the jett iss less tfaanr.TT, . Inr.this scase , 
the buoyant force acts: rrr the negative Yrdirectionv. . In order 
to examine, the effect of' buoyancy, the: cases sinr which Gr/Re 
is non— zero are compared', toe the: case in whichr.Gf/Re :=- 0 . 

For the transverse jet,, the — 0 streamlinesindicates 
the upper boundary af r the flow introduced: by the: jet. No 
flow ever crosses this boundary.. Because- theef low introduced 
by the jet occupies the space: below this: boundary, the flow 
originally in. the- channel, is: flowing through;. aa. smaller area, 
and is therefore accelerated: down stream off the c j e t r . 

A positive value of' Gr/Re: causes the: transverse jet to 
penetrate further into the stream than the jet -for : which 
Gr/Re =• Q, particularly near the jet location.. As a result 
of this larger "obstruction" to the flow in the channel, local 
regions of decelerated and accelerated flow are created near 
the jet (Fig. 17).. A. negative value of Gr/Re restricts the 
penetration of the jet:.. Iir this case, theef low. originally 
in the channel is able: to: surpass the. jet'. inr.a: smooth manner 
(Fig. IffX. 
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Fear ttree Locrpg-.i~ta.idi nail jet* , theeef-f eet .of the buoyant 
force car. the: ffLow.- is.: not: ass pronounced ; as .in the transverse 
jjet.. This? iisE because the - initiallvelocrty of the jet is in 
the- dawnetreranr direction, andd therefore-; the heat added by 
the jet is crmvected' downstream: more -quickly than for the 
transverse jet., Hence , the: local- values of 0 are not as 
great irr. the case of- the" longitudinal 1 jet . Since the 
buoyant force iss proportionail toe 0 , thesbuoyant forces near 
the jet see not- as? great in the-:- case : of : a longitudinal jet. 

The effect of the buoyant: force-' can: be seen by comparing 
the streamlines fox the case: in: which Gr/Re is greater than 
zero to the case in which Gr/Re:— (K. The positive buoyant 
farce causes the streamlines- to - be ■ bent upward with respect 
to the case in. which Gr/Re — 0' (Fig; . 22 } . This is due to 
the velocity in. the positive' Yrdirection caused by this 
buoyant force.. A', negative value-' of :Gf /Re has the opposite 
effect (Fig-.. 232),, causing the: streamlines to be bent 
downward.. 

The effect of' buoyancy on the temperature field for 
both jets can be seen by examining the isotherms presented 
in Figures 19—21 and 25-27. At downstream locations near 
the jet, temperatures are greater. for a positive buoyant 
force than for the case in which there is no buoyant force 
present.. This is- due to the fact that a positive buoyant 
force contributes' to the heat: transfer sin the vertical 
direction' by: increasing the:- vertical 1 velocity component. 
Hence,, the heat added by the jet is dispersed more quickly 
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in. tire vertical! direction, . and' temperatures near the jet 
are higher:.. This results^ inr. anr. increased surface temperature. 
Because afX heat: added' byth'e; jet-must' leave the channel 
through the surface:, . and' thee amounizof heat transferred at 
the surface is proportional! toe thee- surf ace temperature, more 
heat is transferred out off thee channel near the jet for a 
positive buoyant force than. forrae zero buoyant force. 
Therefore the temperaturese ate downstream locations are lower 
for the positive buoyant force;. This effect can be seen by 
examining the surface temperature: data presented in Figure 
20.. This effect: occurs for both:. the longitudinal and trans- 
verse jets,, but is r more pronounced for the transverse jet. 

R. negative- buoyant force has the' opposite effect, causing 
the temperatures near the jet toebe smaller, and the down- 
stream temperatures higher. 
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W.. DISCUSSION 



A. CONVERGENCE: OF KIN EAR- SOLUTION. ‘ FDR? STREAM FUNCTION 

The. validity a£ thee raa-trix of: linearr equations for the 
stream function can be: verified by constructing a conver- 
gence plot tcx determines the value - to: which: the numerical 
solution: wouid’ converge:- iff an inf inite: number of nodes were 
used. This value may then be compared: to: an. exact solution. 

A convergence plott is constructed: by- -plotting the value 
of the numerical soliitian obtained at: a: given point in the 
field of the problem versus 1/N , where: N is the number of 
divisions per side af: the field. By making the number of 
divisions greater,, different values: off the: numerical solution 
are obtained at the same point in the -field. The plot of 
the value of the. solution, versus 1/N. 2 - may- then be extrapo- 
lated to determine, the point- at which, the: plot .crosses the 
1/N 2 = 0 axis.. This value of the solution is then taken as 
the value to which, the numerical solution . would converge if 
an infinite number of divisions were possible. 

In order to apply this method to the linear (zeroth 
order) solution, for the stream function, a" fixed field must 
be chosen. The field chosen was a section of the channel 
of length 5D.. The case considered was a. longitudinal jet 
of width 2H located at the upstream boundary and in the 
center of the channel!. By - utilizing a: grid ~ in which the 
spacing in the X'-direction is five - times: the spacing in. the 
Y-direction, a grid with an equal number of divisions on 
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each edge erf the fredrt is: possibles . Convergence plots were 
ccastrucited at thee points' X ' = 2 '. 5^ , YY== 0^25 and X = 2.5, 

T = G.75.. The mesh, sizes that: wereeusedewere N = 4, N = 8, 
and & = 1 Z- The N' = 44 gxid' was.- too: coarse to yield useful 
ccrarergerce data., Since a grid:, corresponding to N = 16 
requires excessive computer. storage', , only the solutions 
aarresponding to N: = 8! and" N' — 12r. weree-used in the extrapo- 
lation process'.. The convergence, plots: for these two points 
are shown in. Eigures: 311 and 321. . 

An exact solutian at these' two:, points is obtained by 
noting that the surface and bottom: boundary conditions for 
the velocity are the same as those: forr Couette flow between 
flat plates.. The specification of:' a: known flew rate for 
incompressible flaw is equivalent to: specifying a driving 
pressure gradient in the downstream direction. Hence, for 
flow in the channel, far enough downstream so that the effect 
of the jet on the flow has been damped, out, the exact 
solution is the solution for Couette' flow between flat 
plates with an imposed pressure gradient. This solution is 
given in Schlicting [14]. It should be noted, that this 
velocity profile is parabolic in y, and that U (0) =0 and 
U (D) = U s . For an incompressible flow, the pressure 

gradient 1Z. may be determined from a. specified flow rate. 

9 X 

At the downstream boundary, several. assumptions are made. 
It is assumed' that the ef fect of: the: jet .may be neglected, „ 
with the exception of the increased ' flow added by the jet, 
and that v = 0. Tt is also assumed that the velocity profile 
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Eigure 31 — Convergence plot for stream function at 
X = 2.5, Y = 0.75. 
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Figure 22 — Convergence plot for stream function 
at X = 0.25, Y = 2.5. 
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iLss parnboT: irr. in y„ and: that; thee flbwv. rateeis sknown . In 

addit tiicmr.,, the vedrrrriiy-' air y/ ~ O', iss equal itac zero and that 
trite vteHocity.- at y = EC’ iss equal! to: U s , . Sinceethe velocity 
profile- ite parabolic: irr. y;, and' since: theethree conditions 
which. determine- thee coefficients off the: parabola are 
Meatical to the cxmdihions' for Couette- flow with an imposed 
pressure gradient,, tire: velocity prof ile'atithe downstream 
boundary is identical! to: the exact: solution; 

rt was: assumed’. that: the exact- solutionrwas valid at 
X: = 215.. Thee values obtained from:- theeconvergence plot may 
then be compared ter those:- obtained from the exact solution, 
at the point, x: = 2!5f, Y. == 0.25, the exact: solution is 
i|»* =• 0 ..12123 ;• the salution obtained, from- the convergence 
piot is ij>* = a.lZJH.. at: the point" X‘ == 2 : 5 Y = 0.75, the 
exact salution is- \p* =■ O'. 83170; the- solution obtained from 
the convergence plcct is: ip*' — 0.83163; . Asearresult of the 
agreement between, these: values , i.timay- be concluded that 
the linear salution fox which is the zeroth order solu- 
tion for the iterative scheme, is valid. 

B. NUMERICAL STABILITY OF ITERATIVE SCHEME 

By formulating the problem in terms of the linear 
solutions as shown in equations (20) and (21) , the nonlinear 
terms in the energy equation are multiplied by the Peclet 
number Ee,, where Ee — Re Pte. The nonlinear terms in the 
momentum equation: are- multiplied by- Re-; . However, using the 
iterative- procedure described earlier; theenumerical solution 
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■ffcmrr the ener gy equation becomes: unstable at approximately 
Fe = 2.. Since the energy • and! momentum- equations are coupled, 
the e n t-i-re solution. becomes.- unstable: at • Pe = 2. 

Several methods were employed* inr. an - effort to eliminate 
this instability.. The method:, initially; used to calculate 
the. nonlinear.' terms was to: use: thee operators shown in 
Figures' 5,. 5„ and' 7.. Using: this: method , the energy equation 
became unstable at approximately • Pee =- 1. When uncoupled 
from the energy equation, the: momentum equation became 
unstable act Re = 15. 

The method.' at calculating, the nonlinear terms in terms 
of the Jacobian, as described* earlier: was then attempted. 

Using this method, the Peclet: number, at which instability 
occurred in. the energy equation was: raised from 1 to 2 . 

The Reynolds number at which instability occurred in the 
momentum' equation,, uncoupled! from* the energy equation, was 
raised from 15 to 50 . 

Decreasing the grid spacing- was- also attempted. This 
had no noticeable effect on the stability of the system of 
equations ., 

An attempt was made to begin with the stable solution 
for Pe = 2.0, and increase the Peclet number in steps of 
1/25. Using this procedure, instability occurred at Pe = 2.15. 

C. CONCLUSIONS 

Because' all. solutions presented: in this study are based 
cm iterations from the linear; , uncoupled solutions of the 
derived governing equations, the validity of the numerical 
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used: toe obtain these- solutions is of utmost impor- 



tance... Hies vaiidity off th'es operators 'and method of solution 
for: tires lirreair energy equation: has s been demonstrated by its 
frequent usee with successs inr. th'eeliterature [12, 15], The 
results; afs thee convergence' plot: indicate that the linear 
solution of the: Linear momentum* equation is also valid. 

IBased cur theses observations* , if. is: felt that the solutions 
presented^ ih: this: report- aree calculated from an essentially 
sound basis;.. 

Hie effect of buoyancy on:: the: flow patterns and tempera- 
ture distributions has been; demonstrated for the cases 
conrsidereGd'.. However, the* original'. motivation for this 
effort was to: demonstrate the effect .of buoyancy on the flow 
im an. actual', physical channel!. Forrsuch channels, the 
actual values of Reynolds and ’ P'eclet' numbers experienced 
sore, greater; tharr 10^. Because: of:' the restraints on the 
range of Re and'. P.e. due to nonlinear: numerical instability, 
it was nut possible to generate: solutions which are appli- 
cable to actual channels. 

The maximum value of the Peclet number for which results 
were obtained was approximately 2.0. Using the value 
Pr = 6 as a typical value of the Prandtl number for water, 
the maximum Reynolds number for which results for the 
coupled temperature and flow fields "were obtained was 
approximately 0C.3L. 

Rs a. result: of these limitations > some of the assump- 
tions made for the boundary conditions for the longitudinal 



75 ' 



jet may be question eh. . Err. the formulation: of :: the boundary 
condition far: the temperature' for: thee longitudinal jet, it 
was assumed that: mx heat was: conducted: upstream - to the 
X-location af the jefc.. Th±s assumption - , is s valid for 
Pe >> 1. Since the maximum Peclet nnmberr used: in the cou- 
pled problem was Ee = 31. .8E,, this assumption - is £ questionable. 

For the velocity boundary condition:. forrthe longitudinal 
jet, it was- assumed:; that: the distrubancee toe the : flow field 
caused by the get does not: propagate- upstream: to the 
X-lacation of the jet.. This assumption. - , issvalid for 
Re >> L. Because the maximum Reynolds: number:used in the 
solution of the coupled, energy and momentum" equations was 
0.3, this assumption is no longer valid. 

D. RECOMMENDATION S 

Using a finite difference method: off solution, one 
algebraic equation, ie generated for. each.. node ; . In order to 
solve such a system of linear algebraic' equations , two 
basic methods exist.. The first method is: the, direct method, 
which is to solve the simultaneous equations directly , 
generating an exact solution to the system of equations. 

The second method is the indirect method, which arrives at 
a solution by successive improvements of an approximate 
solution. This process' is repeated until' a: desired conver- 
gence criteria is satisfied. 

The basic method of solution used:' in- this study was the 
direct method.. Using this - method/ the- number: of calcula- 
tions required to achieve a solution is known, which is not 



76 



true cf the indirect: method;. The direct: method* is also well 



adapted ta machine csrl'crlartions ♦ Ihr. addition:, the Gaussian 
eliml natron: technique of solution . byy thee direct: method is 
the most efficient method: of' solving- a?, s ye tern." of simultane- 
ous linear algebraic: equations . Using: this: method , however, 
extensive computer storage is required:'.. For: example , a 
finite difference solution for a grid: involving 10 nodes in 
each direction generates: 100 simultaneous: equations for 
each unknown... In. order- tec solve this: problem: by the direct 
method,, a IQCf by- 1CT.GI matrix: must be: stored ‘ for: each unknown. 

On. the basis of this work, the solution of- nonlinear 
algebraic equations using- the iterative: method proposed 
earlier is not feasible when the valuer- of : the nonlinear 
terms becomes large.. For the solution: of-' such a problem, 
the indirect method' of’ solution is ■ recommended. 

One passible method af solving- such' a: system of nonlinear 
algebraic equations is to employ a numerical" relaxation 
scheme. This method is discussed in [12]!. One advantage 
of this method is that once the relaxation pattern is 
established, the need for storing matrices is eliminated. 

Another possible method of solution is to recast the 
problem in the form of vorticity transport 
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where a is the vorticity.. Thet finite: difference-:. forms: of 
these equations cam he solved: for three- value: off theeunknown 
at a node in terms af the values off otherr unknowns; at: 
surrounding nodes. For: instance:,, equation:. (22J; can:. bee 
solved for w in terms of at: surrounding: points; . Thee 
values obtained from the previous iteration:, can:, be: used 
to calculate new values.. By- charrsing: somee appropriate- 
starting point* am iterative. solution. scheme: can:. bee 
established. 

A further recommendation, is that: aa graded: network be 
used in the numerical solution.. This: graded: network should 
have finer divisions near the jet location.. This: would : 
improve the accuracy of the method of' modeling the jet, 
and also provide more grid' paints near, the' jet^ . where the 
maximum changes in flow and' temperature- are experienced. . 
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18 = 17+1 

19 = 17+2 

RETURN 

END 
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